Author: Jose Carlos Molano de Oro
University: Pontificia Universidad Javeriana
Course: Linear Regression Analysis
Semester: 2022-3
Professor: Mario Gregorio Saavedra Rodriguez
Author Email: jose_molano@javeriana.edu.co
Professor Email: saavedrarmg@javeriana.edu.co
El paquete palmerpenguins contiene los datasets:
Se hara el análisis de regresión lineal en base al dataset penguins que posee información de:
se tienen de igual manera para el grupo anterior contiene las medidas de
library(tidyverse)
library(palmerpenguins)
library(reactable)
library(reactablefmtr)
library (DT)
library(Hmisc)
library(GGally)
library(corrplot)
library(ggfortify)
reactable(penguins,rownames = TRUE)
summary(penguins)
## species island bill_length_mm bill_depth_mm
## Adelie :152 Biscoe :168 Min. :32.10 Min. :13.10
## Chinstrap: 68 Dream :124 1st Qu.:39.23 1st Qu.:15.60
## Gentoo :124 Torgersen: 52 Median :44.45 Median :17.30
## Mean :43.92 Mean :17.15
## 3rd Qu.:48.50 3rd Qu.:18.70
## Max. :59.60 Max. :21.50
## NA's :2 NA's :2
## flipper_length_mm body_mass_g sex year
## Min. :172.0 Min. :2700 female:165 Min. :2007
## 1st Qu.:190.0 1st Qu.:3550 male :168 1st Qu.:2007
## Median :197.0 Median :4050 NA's : 11 Median :2008
## Mean :200.9 Mean :4202 Mean :2008
## 3rd Qu.:213.0 3rd Qu.:4750 3rd Qu.:2009
## Max. :231.0 Max. :6300 Max. :2009
## NA's :2 NA's :2
Dado que en las columnas numéricas se tienen valores nulos (NA), estos valores se añadiran con el valor central (mediana) de cada atributo. De igual manera se expresara los años en variable de tipo factor, con fin de conocer cuantos pinguinos se tuvieron en cuenta en el estudio por año
penguins[,3:6]<-penguins[,3:6] %>% mutate_all(~ifelse(is.na(.x), median(.x, na.rm = TRUE), .x))
penguins$year<-as.factor(unlist(penguins$year))
summary(penguins)
## species island bill_length_mm bill_depth_mm
## Adelie :152 Biscoe :168 Min. :32.10 Min. :13.10
## Chinstrap: 68 Dream :124 1st Qu.:39.27 1st Qu.:15.60
## Gentoo :124 Torgersen: 52 Median :44.45 Median :17.30
## Mean :43.92 Mean :17.15
## 3rd Qu.:48.50 3rd Qu.:18.70
## Max. :59.60 Max. :21.50
## flipper_length_mm body_mass_g sex year
## Min. :172.0 Min. :2700 female:165 2007:110
## 1st Qu.:190.0 1st Qu.:3550 male :168 2008:114
## Median :197.0 Median :4050 NA's : 11 2009:120
## Mean :200.9 Mean :4201
## 3rd Qu.:213.0 3rd Qu.:4750
## Max. :231.0 Max. :6300
datatable(as.matrix(sapply(penguins[,3:6],function(x) mean(x, na.rm=TRUE))))
Las medidas promedio en toda la población de la muestra de pinguinos en (mm) y su masa en (g) son:
c=cor(penguins[,3:6])
y=as.data.frame(c)
y[y==1]<-" "
y <- mutate_all(y, function(x) as.numeric(as.character(x)))
reactable(as.data.frame.array((y)),
defaultColDef = colDef(
style = highlight_min_max(as.data.frame.array((y)))))
corrplot(cor(penguins[,3:6],use = "complete.obs"),method="number")
corrplot(cor(penguins[,3:6],use = "complete.obs"),method="circle")
Se puede apreciar que los atributos que están más correlacionados positivamente entre si por orden son:
Para la construcción del modelo de regresión se construira un modelo entre las variables de masa del cuerpo (body_mass_g) vs. longitud de aleta (flipper_length_mm) y longitud de aleta (flipper_length_mm) y altura del pico (bill_depth_mm).
Para las variables masa del cuerpo (body_mass_g) vs. longitud de pico (bill_length_mm) se descarta la construcción del modelo dado que su estadistico de correlación esta muy cercano al 0.5, es decir que no se podria concluir bien si sus valores están directamente correlacionados o no.
ggscatmat(penguins, columns=3:6, color="species")
A partir del gráfico anterior se puede apreciar lo siguiente:
El gráfico de dispersión de las variables longitud de aleta (flipper_length_mm) vs masa del cuerpo (body_mass_g) tiene un comportamiento muy similar a una función lineal con pendiente positiva teniendo en cuenta todas las especies. A mayor peso del pinguino mayor sera su longitud de aleta.
El gráfico de dispersión de las variables longitud de pico (bill_length_mm) vs masa del cuerpo (body_mass_g) tiene un comportamiento muy similar a una función lineal con pendiente positiva teniendo en cuenta solo las especies de tipo Adelie (0.67) y Gentoo (0.59). A mayor peso del pinguino mayor sera su longitud de pico.
El gráfico de dispersión de las variables longitud de pico (bill_length_mm) vs longitud de aleta (flipper_length_mm) tiene un comportamiento muy similar a una función lineal con pendiente positiva teniendo en cuenta solo las especies de tipo Gentoo (0.7). A mayor longitud de pico mayor sera la longitud de aleta del pinguino.
Se proceden a construir algunos modelos de regresión lineal teniendo en cuenta las variables que tienen mas correlación con los gráficos anteriores.
ggplot(penguins, aes(y = flipper_length_mm,
x = body_mass_g,
)) +
geom_point() +
stat_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
\[ \textrm{flipper_length_mm}=\beta_0+\beta_1\times \text{body_mass_g}+e \]
\[ e_i\sim\mathcal{N}(0,\sigma^2) \] \[ f(x)=m\cdot x+b+e \]
model1<-lm(flipper_length_mm~body_mass_g,data=penguins)
model1
##
## Call:
## lm(formula = flipper_length_mm ~ body_mass_g, data = penguins)
##
## Coefficients:
## (Intercept) body_mass_g
## 136.71105 0.01528
\[ \widehat{\textrm{flipper_length_mm}}_i=\widehat\beta_0+\widehat\beta_1\times \text{body_mass_g}_i+e_i \] \[ =1.367\times 10^2+1.528\times10^{-2}\times \text{body_mass_g}_i+e_i \] \[ e\sim\mathcal{N}(0,\sigma^2) \] \[ \widehat{\textrm{flipper_length_mm}}_i=\widehat\beta_0+\widehat\beta_1\times \text{body_mass_g}_i \] \[ =1.367\times 10^2+1.528\times10^{-2}\times \text{body_mass_g}_i \] \[ e_i=\textrm{flipper_length_mm}_i-\widehat{\textrm{flipper_length_mm}}_i \]
summary(model1)
##
## Call:
## lm(formula = flipper_length_mm ~ body_mass_g, data = penguins)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.7543 -4.8130 0.9608 5.0946 16.6487
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.367e+02 1.990e+00 68.68 <2e-16 ***
## body_mass_g 1.528e-02 4.655e-04 32.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.894 on 342 degrees of freedom
## Multiple R-squared: 0.759, Adjusted R-squared: 0.7583
## F-statistic: 1077 on 1 and 342 DF, p-value: < 2.2e-16
A partir de lo anterior se puede observar:
ggplot(penguins, aes(x = body_mass_g,
y = flipper_length_mm,
)) +
geom_point() +
stat_smooth(method=lm, se = TRUE, color = "red")
## `geom_smooth()` using formula 'y ~ x'
confint(model1)
## 2.5 % 97.5 %
## (Intercept) 132.79589823 140.6261928
## body_mass_g 0.01436252 0.0161937
penguins$fitted_values<-model1$fitted.values
penguins$residuals<-model1$residuals
ggplot(penguins, aes(x = fitted_values, y = residuals)) +
geom_point(aes(color = residuals)) +
scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
geom_hline(yintercept = 0) +
geom_segment(aes(xend = fitted_values, yend = 0), alpha = 0.2) +
theme_bw() +
theme( legend.position = "none")
Los residuos se distribuyen de manera aleatoria entorno a la recta \(y=0\), por lo cual se acepta la linealidad.
ggplot(penguins, aes(x = residuals)) +
geom_histogram(aes(y = ..density..),bins=90,fill='white',color=4) +
theme_light()
autoplot(model1)
shapiro.test(model1$residuals)
##
## Shapiro-Wilk normality test
##
## data: model1$residuals
## W = 0.98335, p-value = 0.0005215
Dado que en mi test anterior, mi \(p-\text{value}<0.05\), se rechaza la hipótesis nula por lo cual los residuos no presentan un comportamiento normal.
ggplot(penguins, aes(x = fitted_values, y = residuals)) +
geom_point(aes(color = residuals)) +
scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
geom_segment(aes(xend = fitted_values, yend = 0), alpha = 0.2) +
geom_smooth(se = FALSE, color = "firebrick") +
geom_hline(yintercept = 0) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
AdelieGentoo<-subset(penguins,species!="Chinstrap")
ggplot(AdelieGentoo, aes(y = bill_length_mm,
x = body_mass_g,
)) +
geom_point() +
stat_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
\[ \textrm{bill_length_mm}=\beta_0+\beta_1\times \text{body_mass_g}+e \]
\[ e\sim\mathcal{N}(0,\sigma^2) \] \[ f(x)=m\cdot x+b+e \]
model2<-lm(bill_length_mm~body_mass_g,data=AdelieGentoo)
model2
##
## Call:
## lm(formula = bill_length_mm ~ body_mass_g, data = AdelieGentoo)
##
## Coefficients:
## (Intercept) body_mass_g
## 19.230360 0.005441
\[ \widehat{\textrm{bill_length_mm}}_i=\widehat\beta_0+\widehat\beta_1\times \text{body_mass_g}_i+e_i \] \[ =1.923+5.441\times10^{-3}\times\text{body_mass_g}+e_i \] \[ e_i\sim\mathcal{N}(0,\sigma^2) \] \[ \widehat{\textrm{bill_length_mm}}_i=\widehat\beta_0+\widehat\beta_1\times \text{body_mass_g}_i \] \[ =1.923+5.441\times10^{-3}\times\text{body_mass_g} \]
\[ e_i=\textrm{bill_length_mm}_i-\widehat{\textrm{bill_length_mm}}_i \]
summary(model2)
##
## Call:
## lm(formula = bill_length_mm ~ body_mass_g, data = AdelieGentoo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.5720 -1.7594 -0.0742 1.6852 7.4499
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.923e+01 7.977e-01 24.11 <2e-16 ***
## body_mass_g 5.441e-03 1.815e-04 29.98 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.508 on 274 degrees of freedom
## Multiple R-squared: 0.7664, Adjusted R-squared: 0.7655
## F-statistic: 898.9 on 1 and 274 DF, p-value: < 2.2e-16
ggplot(AdelieGentoo, aes(x = body_mass_g,
y = bill_length_mm,
)) +
geom_point() +
stat_smooth(method=lm, se = TRUE, color = "red")
## `geom_smooth()` using formula 'y ~ x'
confint(model2)
## 2.5 % 97.5 %
## (Intercept) 17.659872328 20.800848515
## body_mass_g 0.005083984 0.005798569
AdelieGentoo$fitted_values<-model2$fitted.values
AdelieGentoo$residuals<-model2$residuals
ggplot(AdelieGentoo, aes(x = fitted_values, y = residuals)) +
geom_point(aes(color = residuals)) +
scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
geom_hline(yintercept = 0) +
geom_segment(aes(xend = fitted_values, yend = 0), alpha = 0.2) +
theme_bw() +
theme( legend.position = "none")
Los residuos se distribuyen de manera aleatoria entorno a la recta \(y=0\), por lo cual se acepta la linealidad.
ggplot(AdelieGentoo, aes(x = residuals)) +
geom_histogram(aes(y = ..density..),bins=90,fill='white',color=4) +
theme_light()
autoplot(model2)
shapiro.test(model2$residuals)
##
## Shapiro-Wilk normality test
##
## data: model2$residuals
## W = 0.99581, p-value = 0.6686
Dado que en mi test anterior, mi \(p-\text{value}>0.05\), se rechaza la hipótesis alternativa por lo cual los residuos presentan un comportamiento normal.
ggplot(AdelieGentoo, aes(x = fitted_values, y = residuals)) +
geom_point(aes(color = residuals)) +
scale_color_gradient2(low = "blue3", mid = "grey", high = "red") +
geom_segment(aes(xend = fitted_values, yend = 0), alpha = 0.2) +
geom_smooth(se = FALSE, color = "firebrick") +
geom_hline(yintercept = 0) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'